home *** CD-ROM | disk | FTP | other *** search
/ CD ROM Paradise Collection 4 / CD ROM Paradise Collection 4 1995 Nov.iso / program / nrpas13.zip / MEDFIT.PAS < prev    next >
Pascal/Delphi Source File  |  1991-04-29  |  2KB  |  65 lines
  1. PROCEDURE medfit(x,y: glndata; ndata: integer; VAR a,b,abdev: real);
  2. (* Programs using routine MEDFIT must define the type
  3. TYPE
  4.    glndata = ARRAY [1..ndata] OF real;
  5. in the main routine. MEDFIT also assumes that the instructions
  6. at the beginning of ROFUNC have been carried out so that arrays
  7. x,y and arr, and real variables aa and abdevt exist as globally
  8. defined variables. *)
  9. LABEL 3;
  10. VAR
  11.    j: integer;
  12.    sy,sxy,sxx,sx,sigb,f2,f1,f: real;
  13.    del,chisq,bb,b2,b1: real;
  14. BEGIN
  15.    sx := 0.0;
  16.    sy := 0.0;
  17.    sxy := 0.0;
  18.    sxx := 0.0;
  19.    FOR j := 1 TO ndata DO BEGIN
  20.       sx := sx+x[j];
  21.       sy := sy+y[j];
  22.       sxy := sxy+x[j]*y[j];
  23.       sxx := sxx+sqr(x[j])
  24.    END;
  25.    del := ndata*sxx-sx*sx;
  26.    aa := (sxx*sy-sx*sxy)/del;
  27.    bb := (ndata*sxy-sx*sy)/del;
  28.    chisq := 0.0;
  29.    FOR j := 1 TO ndata DO BEGIN
  30.       chisq := chisq+sqr(y[j]-(aa+bb*x[j]))
  31.    END;
  32.    sigb := sqrt(chisq/del);
  33.    b1 := bb;
  34.    f1 := rofunc(b1);
  35.    IF (f1 > 0.0) THEN BEGIN
  36.       b2 := bb+abs(3.0*sigb)
  37.    END ELSE BEGIN
  38.       b2 := bb-abs(3.0*sigb)
  39.    END;
  40.    f2 := rofunc(b2);
  41.    WHILE ((f1*f2) > 0.0) DO BEGIN
  42.       bb := 2.0*b2-b1;
  43.       b1 := b2;
  44.       f1 := f2;
  45.       b2 := bb;
  46.       f2 := rofunc(b2)
  47.    END;
  48.    sigb := 0.01*sigb;
  49.    WHILE (abs(b2-b1) > sigb) DO BEGIN
  50.       bb := 0.5*(b1+b2);
  51.       IF ((bb = b1) OR (bb = b2)) THEN GOTO 3;
  52.       f := rofunc(bb);
  53.       IF ((f*f1) >= 0.0) THEN BEGIN
  54.          f1 := f;
  55.          b1 := bb
  56.       END ELSE BEGIN
  57.          f2 := f;
  58.          b2 := bb
  59.       END
  60.    END;
  61. 3:   a := aa;
  62.    b := bb;
  63.    abdev := abdevt/ndata
  64. END;
  65.